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We study the shapes of human red blood cells using continuum mechanics. In particular, we model 
the crenated, echinocytic shapes and show how they may arise from a competition between the 
bending energy of the plasma membrane and the stretching/shear elastic energies of the membrane 
skeleton. In contrast to earlier work, we calculate spicule shapes exactly by solving the equations of 
continuum mechanics subject to appropriate boundary conditions. A simple scaling analysis of this 
competition reveals an elastic length Aej which sets the length scale for the spicules and is, thus, 
related to the number of spicules experimentally observed on the fully developed echinocyte. 

I. INTRODUCTION 

The normal resting shape of the human red blood cell (RBC or erythrocyte) is a flattened biconcave disc (discocyte) 
about 8 ^m in diameter. Treatment of erythrocytes in vitro with a variety of amphipathic agents is known to 
transform this shape systematically and reversibly into various other shapes such as echinocytes (crenated shapes) 
and stomatocytes (cup-like shapes), which are further subdivided into classes labelled I, II, III (Brecher and Bessis, 
1972; Bessis, 1973; Chailley et al. 1973; Mohandas and Feo, 1975). In particular, the echinocyte III is a more-or-less 
spherical body covered evenly with 25-50 rounded protuberances, which we shall call spicules. 

These shape transformations — from the discocyte to the stomatocyte on one side and from the discocyte to the 
echinocyte on the other — have been studied experimentally for more than 50 years. Understanding these shapes 
and shape transformations is a classic problem in cell biology; over the past three decades it has also attracted the 
attention of physicists. What makes this problem so intriguing is that the structure of the red cell is remarkably 
simple. To a good approximation it is simply a bag of fluid, a concentrated solution of hemoglobin surrounded by a 
thin macroscopically homogeneous membrane. Thus, the diverse resting shapes which it exhibits can be thought of as 
(locally) stable mechanical equilibrium shapes of the membrane. This membrane is composite (Steck, 1989; Alberts 
et al., 1994). On the outside is the plasma membrane, a self-assembled fluid bilayer with a thickness of about 4 nm, 
composed of a complex mix of phospholipids, cholesterol, and dissolved proteins. Lipids such as phosphatidylcholine, 
sphingomyelin, and glycophospholipids, which are neutral at physiological pH, are common in the outer leaflet, while 
phosphatidylserine and phosphatidylethanolamine predominate in the inner leaflet (Gennis, 1989; Alberts et al., 
1994). Due to the negative charge of phosphatidylserine, located in the inner leaflet, there is a significant difference 
in charge between the two leaflets of the bilayer. Inside the plasma membrane but linked to it by protein anchors is a 
thin, two-dimensionally cross-linked protein cytoskeleton, the membrane skeleton (Steck, 1989; Bennett, 1990). The 
skeleton is an hexagonally linked net, each unit of which is a fllamentous spectrin tetramer with an extended length 
of about 200 nm. The spectrin tetramer is negatively charged at physiological pH. Thermal fluctuations reduce the 
projected distance between the vertices of the net to about 76 nm, while the offset between the plasma membrane and 
the spectrin network is 30-50 nm (Byers and Branton, 1985; Liu et al., 1987). Functionally, the plasma membrane 
serves as an osmotic barrier, passing water with relative ease but controlling, via a system of pumps and channels, 
the flow of ions and larger solute molecules. The membrane skeleton has the role of supporting and toughening the 
plasma membrane, which would otherwise disintegrate in circulatory shear flow. 

While the RBC membrane is certainly heterogeneous at the length scale of individual lipid molecules, lipid patches 
or the basic spectrin tetramer, on scales longer than 100 nm it is reasonably homogeneous in its properties, and it 
may make sense to treat it as a mechanical continuum. In this spirit and following early work by Canham (1970), 
Hclfrich (Helfrich, 1973; Deuling and Hclfrich, 1976), and others, we imagine replacing the (fluid-phase) plasma 
membrane by an ideal uniform two-dimensional surface characterized by a bending rigidity taken for the RBC to 
be (Waugh and Bauserman, 1995; Strey et al., 1995) 2.0 x 10~^^ J (roughly 50 fesTroom, where ks is Boltzmann's 
constant), a spontaneous curvature Co, which recognizes the outside-inside asymmetry of the leaflet composition, and 
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an arca-difFerence-elasticity (ADE) term (Helfrich, 1973; Evans, 1974; Svctina et al., 1985; Bozic et al., 1992; Wicse 
et al., 1992; Miao et al., 1994), which reflects the fact that a difference in relaxed area, AAq, between the outer and 
inner leaflets produces a "bilayer couple" (Sheetz and Singer, 1974) tending to bend the membrane away from the 
larger-area side (like a bimetallic strip). We assign to the modulus R associated with the ADE term a value of (Waugh 
and Bauserman, 1995) 1.27 x 10~^^ J. These terms taken together constitute what is called the ADE model of the 
plasma membrane. In addition, we represent the RBC membrane skeleton as a two-dimensional elastic continuum 
characterized by a bulk stretching modulus Ka and a shear modulus fx with w 2.5 x 10~^ J/m^ and Ka « 3/x 
(see discussion in Sec. V). These two mechanical subsystems, constrained to the same mathematical surface, will 
constitute our model of the RBC shape at length scales above 100 nm (see Sec. II for full details). 

The basic hypothesis of the mechanical approach is that the observed RBC shapes must be shapes which minimize 
the energy subject to appropriate constraints on volume and area, i.e., that all the observed RBC shapes must emerge 
as local energy minima of this model and that the observed shape transformations must come about as a response of 
the minimum-energy shapes to changes in the "control parameters" which characterize the model. There are relatively 
few such control parameters: In addition to the known area and volume of the RBC, we include in the list the three 
mechanical moduli noted above, the spontaneous curvature Co, the relaxed area difference AAq between the bilayer 
leaflets, and parameters describing the effective relaxed size and shape of the membrane skeleton. If all these control 
parameters were; known or easily measurable, we could predict the stable RBC shapes using continuum mechanics. 
The fact that they are not restricts us, as we shall see, to somewhat more generic predictions. 

This program has already been carried out, albeit in somewhat restricted form, with respect to the discocyte 
and stomatocyte shapes. The shapes and shape transitions of laboratory-prepared fluid-phase phospholipid vesicles 
(lacking any skeleton) have been successfully described by the ADE model (Miao et al., 1994; Dobereiner et al., 
1997). In particular, plausible discocytic and stomatocytic shapes do occur as energy-minimizing shapes of the ADE 
model. Furthermore, the bilayer-coupic mechanism (Sheetz and Singer, 1974) accounts qualitatively for many of 
the chemically-induced discocyte-stomatocytc transitions observed in the laboratory, in that stomatocytogenic agents 
tend to segregate to the inner leaf of the bilayer, thus making AAq negative and favoring the invaginated stomatocyte 
shape (see more below). There has been discussion in the literature of the way in which introducing the elastic 
energy of the membrane skeleton modifies or improves the detailed prediction for the shape of the normal resting 
discocyte (Zarda et al., 1977; Evans and Skalak, 1980). However, most treatments in the literature tend to regard 
the cytoskeletal elastic energy as providing only perturbative modification of the basic shapes which emerge from the 
pure-bending-energy model. We shall argue that this point of view fails completely in describing echinocyte shapes. 

Until recently, it was a major difficulty for a fully mechanical picture of the RBC shapes that cchinocytic shapes 
have not been found in the catalog of minimum-energy shapes of the ADE model; however, several authors have now 
pointed out what we believe is the correct resolution of this problem (Waugh, 1996; Iglic, 1997; Iglic et al., 1998a, 
b; Wortis, 1998): A pure- lipid vesicle with a sufficiently positive spontaneous curvature (or, equivalently, sufficiently 
large positive AAq) adopts a vesiculated or budded shape, in which one or more quasi-spherical buds appears on the 
outside of the main body of the vesicle, attached to it via a narrow neck or necks (sec Fig. 1). Such necks are allowed 
as low-energy structures for pure-lipid vesicles, because the bending energy depends on the sum of the two principal 
curvatures (see Sec. II), so that large curvatures of opposite sign (characteristic of necks) can still lead to a small mean 
curvature and, thus, to low energy (Fourcade et al., 1994). But — and this is the crucial point small necks involve 
large elastic strains in the membrane skeleton and are, thus, inhibited by elastic energy cost. It is not surprising 
that they are not seen for undamaged RBC's. What our calculations show is that the formation of echinocytes is 
driven by positive Cq (equivalently, positive AAq) and that the spiculated shapes seen in experiments can arise from 
a competition between the bending energy of the plasma membrane, which by itself would promote budding, and the 
stretching and shear elastic energies of the membrane skeleton, which prevent it. To make this more precise, we note 
that the ratio of the bending modulus to the shear modulus (or the stretching modulus) defines an elastic length scale, 



which yields A^i ~ 0.28 /im. We shall sec below that this number sets the scale of spicule size and, thus, fixes the 
spicule density of the fully developed echinocyte (echinocyte III). 

The success of this approach strengthens the hypothesis that RBC shapes can be understood on the basis of 
continuum mechanics; however, in themselves, such calculations cannot answer the question, "Why does a particular 
shape occur under given experimental conditions?" Such a question has two parts. Why do the control parameters 
have the values they do? And, how does this set of control parameters lead to the observed shape? The second 
part is the proper domain of continuum mechanics. The first part is primarily biological or biochemical in content 
and revolves around the mechanisms which control the lipid composition of the leafiets of the plasma membrane, 
local electrostatic effects, the composition and environment of the spectrin network and its coupling to the plasma 
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membrane, and any other determinants of Co, AAq, the elastic moduli, and other control parameters. Hopefully, this 
separation of the question is useful. 

The plan of this paper is as follows. In the remainder of this section, we briefly review some history of echinocyte 
shapes and shape calculations. Section II presents the details of the continuum-mechanics model. Section III discusses 
the theory required to solve the model numerically for echinocyte shapes, including a full account of the scaling 
argument suggested above. Section IV presents our results for the predicted echinocyte shapes and spicule density. 
Section V contains discussion of our results and further speculations. 

The discocyte-to-echinocyte transformation was precisely identified by Ponder (1948, 1955) and has since been 
studied extensively by many other authors, as reviewed, for example, by Bessis (1973), Lange et al. (1982), and 
Steck (1989). Although this transformation can be driven in many ways, the final spiculated shapes produced are 
apparently the same (Smith et al., 1982; Mohandas and Feo, 1975). This quasi-universal behavior is evidence for the 
kind of mechanism suggested above, in which quite diverse biochemical processes may set (a few) important control 
parameters to the same or similar values. In the same spirit, we note that RBC ghosts (formed by hemolysis) behave 
quite similarly to intact red cells (Lange et al., 1982). 

The simplest picture of this type — consistent with some (but not necessarily all) experimental findings — is the 
so-called bilayer-couple hypothesis of Sheetz and Singer (1974), which focuses attention on the control parameter 
A^o- Thus, adding exogenous phospholipids to the exterior solution is observed to promote echinocyte formation 
(Ferrell et al., 1985; Christiansson et al., 1985). This is explained in the bilayer-couple picture by arguing that 
such addition is expected to produce an area increase in the outer leaflet only, since phospholipids do not readily 
flip-flop from one leaflet to the other. This increases the area difference AAq and is expected to promote echinocyte 
formation. Similarly, cholesterol does equilibrate across the bilayer but apparently prefers the outer leaflet, presumably 
because of the particular lipids there present. Thus, adding cholesterol will tend to increase A^o and promote 
echinocytosis, while depicting cholesterol will tend to decrease A^o and to promote stomatocytosis (Lange and 
Slayton, 1982). More generally, anionic amphipaths all tend to produce echinocytcs, while cationic amphipaths tend 
to produce stomatocytes (Deuticke, 1968; Weed and Chailley, 1973; Mohandas and Feo, 1975; Smith et al., 1982). The 
bilayer-couple explanation of these observations is that, when incorporated into the plasma membrane, the cationic 
compounds segregate preferentially to the inner leaflet and the anionic compounds, to the outer leaflet because of 
the predominantly negative charge of the inner- leaflet lipids. Outer-leaflet segregation increases A^o, thus promoting 
echinocytosis; conversely, inner-leaflet segregation promotes stomatocytosis. Time-dependent effects have even been 
seen, where a molecule initially intercalates into the outer leaflet, causing echinocytosis, but then slowly migrates to 
the inner leaflet, following the electrostatic gradient, causing return to the discocyte and subsequent stomatocytosis 
(Isomaa et al., 1987). 

A similar hypothesis is made to explain the observed tendency of ghosts to become echinocytes at high salt con- 
centrations but stomatocytes at low salt concentrations (Lange et al., 1982). The argument relies on electrostatics. 
Cationic species neutralize the negative lipid charges on the inner leaflet and salt decreases the Debye length, more 
effectively screening such charges. Both effects lead to a contraction of the inner leaflet and consequent increase in 
AAq, promoting crenation. 

Other observations can be linked to the bilayer-couple effect but somewhat less directly. Thus, it is found that 
hemolyzed echinocytes produce discocytic ghosts (Lange et al., 1982) and that electroporation suppresses shape 
changes (Schwarz et al., 1999a, b). To explain these effects, it is argued that hemolysis and electroporation both 
involve perforation of the membrane and resulting in contact via the pore surface between the lipids in the inner and 
outer leaflets. This provides a mechanism for lipid transport which may reduce the magnitude of AAq by partial 
symmetrization of the lipid composition of the leaflets, thus (equivalently) reducing Cq. Similarly, it is known that 
ATP depletion drives discocytes to crenate (Nakao et al., 1960, 1961, 1962; Backman, 1986). The hypothesis here is 
that ATP is required in some way to maintain the lipid asymmetry of the leaflets (related to Co), perhaps for the 
operation of ATP-driven translocases (Steck, 1989). 

Some effects are not readily explained by the bilayer-couple mechanism. Weed and Chailley (1972), and Gedde et 
al. (1995, 1997a, 1997b, 1999) report that RBC shape can be controlled experimentally by varying the external pH, 
with high pH promoting echinocyte formation and low pH, stomatocyte formation (the effect of proximity to a glass 
surface in promoting echinocytosis (Furchgott and Ponder, 1940) is probably related to this effect). The mechanism 
for this effect does not seem to be well established. Some authors argue for a mechanism involving membrane proteins. 
Band 3 has been speciflcally implicated (Gimsa and Ried, 1995; Jay, 1996; Gimsa, 1998; Hagerstrand et al., 2000). 
For example, it is known that high pH induces dissociation of ankyrin and Band 3, which plays a role in the linkage 
of the membrane skeleton to the bilayer (Low et al., 1991). Indeed, it was shown that in vitro the membrane skeleton 
expands at high pH and contracts at low pH (Elgsaeter et al., 1986; Stokke et al., 1986). It has recently been 
proposed that pH change may induce conformational transformation in band 3 protein, leading to a change in A^g- 
Similarly, Wong (1994, 1999) has proposed that the folding of the spectrin in the cytoskeletal net is controlled by 
a conformational change of the Band 3 protein. Such mechanisms might shift the emphasis of shape determination 
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from the dominantly lipid-related control parameters, Cq and AAq, to the elastic parameters and the relaxed shape 
of the membrane skeleton. We shall comment briefly on these issues in Sec. V. 

Numerous shape calculations have previously been done using variants of the model we describe in Sec. II. Early 
calculations (Fung and Tong, 1968; Zarda et al., 1977; Evans and Skalak, 1980) omitted the ADE term (k — 0), set 
Co = 0, and concentrated entirely on the shape of the resting discocyte and the modifications caused by osmotic 
swelling. They solved numerically the full mechanical shape equations. Landman (1984) treated a viscoelastic shell 
surrounding a viscous droplet and modeled the formation of protrusions as sudden local addition of shell material. 
Waugh (1996) studied a full ADE model with local shear elasticity (but treating the compressibility modulus as 
infinite) and studied the instability of a flat membrane to formation of a local "bump" (axisymmetric spicule) for 
positive values of Co or AAq or both. He assumed a specific spicule-shape parametrization, focussed on a single 
spicule, calculated its energy, and in this way was able to discuss generically the conditions under which spicule 
formation might become energetically favorable. Most recently, Iglic and collaborators (Iglic, 1997; Iglic et al., 1998a, 
b) have used the same model as Waugh (1996) but parametrize the spicule more simply, using a hemispherical cap 
on a cylindrical body and joining this to a spherical surface via a base shaped like a toroidal section. They correctly 
identify the emergence of echinocyte shapes as a competition between the bending energy and the skeletal elasticity. 
They model the echinocyte as a collection of spicules attached to a spherical body and then determine the radius 
of the sphere, the parameters of the spicules, and the preferred number of spicules by minimizing the mechanical 
energy subject to constraints on cell area and volume. In this way, they are able to estimate the expected number of 
spicules and their shape as a function of the relaxed area difference AAq. 

The present paper builds on the work of Waugh and Iglic: We extend the model to include the area modulus Ka 
of the membrane skeleton. We calculate spicule shape for the first time from the Euler-Lagrange equations (these 
calculations are closely related to the earlier work of Zarda et al. (1977)). Although we continue with Iglic et al. 
(1998a, b) to treat the echinocyte as a sphere decorated with spicules, we deal seriously (although still approximately) 
with the mechanical boundary conditions which must be met where the spicule joins the sphere. 



II. THE MODEL 



The central assumption of our work is that the red blood cell assumes a shape that minimizes its overall membrane 
energy subject to the appropriate constraints. The RBC membrane is a composite of the plasma membrane and the 
cytoskeletal network; correspondingly, we take the membrane energy to be a sum of two terms, 

F^Fh + F^i, (2) 

the bending energy. 

Ft - ^K, J dA(Ci +C2- CoY + |^(AA - AAo)^ (3) 

associated with the bilayer and an elastic energy of stretching and shear, 

F,i = J dAo(AiA2 - 1)2 + 1 1 dAo (^^ + ^ - 2^ , (4) 

associated with the skeleton. In doing this, we treat the plasma membrane as incompressible. An estimate of its 
stretching modulus is j^j^^^^^'^'^^ ^ 10~^ J/m^, comparable to that of a pure-phospholipid bilayer and some four orders 
of magnitude larger than that of the skeleton. Correspondingly, we ignore any bending rigidity associated with the 
isolated membrane skeleton (the dimensional estimate ^ x (thickness) ^ leaves it two orders of magnitude 

smaller than that of the bilayer) . 

Eq. ^ for the plasma membrane's contribution to the overall membrane energy is the now-standard ADE Hamiltonian 
(Svetina et al., 1985; Bozic et al., 1992; Wiese et al., 1992; Miao et al., 1994). The first term was originally proposed 
by Helfrich (Helfrich, 1973; Deuling and Helfrich, 1976). Ci and C2 are the two local principal curvatures, Co is the 
spontaneous curvature, and the integral is over the membrane surface. The two leaflets of a closed bilayer of fixed 
inter leaflet separation D are required by geometry to differ in area by an amount, 

AA^D J dA(Ci +C2). (5) 

In calculations, we shall take D « 3nm, corresponding to the distance between the neutral surfaces of the leaflets. 
The second term in Eq. || is the so-called area-difference-elasticity energy and represents the cost in stretching (or 
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compressional) energy of the individual leaflets necessary to force the change from the relaxed area difference A^o so 
as to conform to this geometric requirement. This effect occurs because a strong hydrophobic barrier prevents lipids 
in one leaflet from "flip-flopping" passively to the other on the time scales of mechanical shape changes. A is the area 
of the plasma membrane. The modulus k associated with the ADE term is generically of the same scale as kj,. Note, 
finally, that we can substitute Eq. || to rewrite Eq. ||, 



Fb = ^Kb 



J dA{Ci + C2Y + ^ (/ dA{Ci + C2) j - 2Co j dA{Ci + C2 



const., (6) 



where 



Co^Co + ^, (7) 

a — R/ Kbi and the constant term is shape independent. This shows that, in determining minimal shapes at given values 
of the control parameters, Co and H^Aq do not enter independently but only in the form of an effective spontaneous 
curvature Cq or, equivalently, an effective relaxed area difference, 

A3o-AAo + ^, (8) 
7ra 

which we shall often quote in the dimensionless, reduced form 

A^o = ^ = Aao + ^. (9) 
A na 

Note for future reference that increasing Aao by 0.01 is equivalent to increasing Co by 6.7/im~^. 

Eq. H measures the elastic-energy cost of the spectrin network. It depends both on the relaxed shape of the 
membrane skeleton and on the way it is actually distributed over the membrane surface. (The notion of a relaxed 
shape is, of course, somewhat nominal, since removing the skeleton from the plasma membrane-which can certainly 
be done (Sheetz, 1979) -would radically change its local biochemical environment, in a way which would modify its 
shape and elastic constants). In this redistribution, each element of the network will be stretched or compressed. The 
quantities Ai and A2 are the local principal extension ratios of each membrane element (Evans and Skalak, 1980). Kq. 
and fi are the (two-dimensional) moduli for stretch and shear, respectively, as introduced in Sec. I, and the integrals 
are over the undeformed shape. In writing Eq. ^, we are assuming, as appears to be typical, that the membrane 
skeleton does not plastically deform in the course of the experimental shape transformation being described. Note 
that Eq. ^ makes a particular choice of terms in the elasticity beyond those quadratic in the weak-deformation 
parameters, = Ai — 1. For the large elastic strains which are present at narrow necks and for small spicules, these 
nonlinearities do play a role. As far as we are aware, RBC elasticity has not been measured well enough to produce a 
clear preference for a particular form of the nonlinearities. Various authors have used Eq. ^ for the elastic energy in 
other contexts with apparent success; however, alternative forms of the elasticity have also been proposed for dealing 
with problems involving large local deformations (Evans and Skalak, 1980; Discher et al., 1994). We shall make some 
additional comments in Sec. V. 

One attractive feature of Eq. ^ is that the effect of changing the size of the relaxed membrane skeleton by a pure, 
uniform dilation is particularly simple. Suppose that the skeleton is decreased in linear scale by a factor b. This 
means that undeformed areas are reduced by a factor 1/6^, while the extension ratios Ai^2 are each increased by a 
factor b. Note that the integrand of the shear term in Eq. ^ is invariant under this change but the factor A1A2 in the 
stretching integrand increases by . Only the quadratic part of the stretching term influences the membrane shape 
because J dAoAiA2 is just the membrane area A, which is fixed by the incompressibility of the plasma membrane. 
It follows that the effect of decreasing the size of the skeleton is completely equivalent to keeping the size fixed but, 
instead, changing the elastic moduli according to 

< = b^K„ (10) 
fi' = b-^fi, (11) 

i.e., this "prestress" in the membrane skeleton is completely equivalent to no prestress, a harder stretching modulus, 
and a softer shear modulus. It is not known with certainty what "prestress" exists in the typical RBC skeleton. 
Computer simulation has suggested that the relaxed skeleton may be as much as 10-20% smaller in area than the 
plasma membrane (Boal, 1994). On the other hand, the experiment by Svoboda et al. (1992) shows that isolated 
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skeletons are expanded. In our calculations, we shall assume zero net prestress (relaxed skeletal area equal to RBC 
area, i.e., 7 = A/Aq = 1). Any deviation from this would result in modified, effective elastic constants according to 
Eqs. |l^ and |ll]. 

In our calculations, wc will identify shapes which arc minima of the (free) energy defined by Eqs. |^-^ subject to 
constraints of fixed membrane area A (bilayer incomprcssibility) and volume V (since the RBC volume is typically 
set by osmotic balance). Constrained minimization is achieved variationally by introducing the functional, 

<i>{Y.,P) = F + J:A- PV, (12) 

where S and P are Lagrange multipliers used to enforce the surface-area and volume constraints. P is the pressure 
difference across the bilayer, while S has the dimension of a surface tension. Making $ stationary with respect to 
variations of membrane shape and cytoskeletal distribution leads to a set of coupled Euler-Lagrange equations. These 
equations can then be solved to give the shapes of mechanical equilibrium. In principle, such shapes are expected to 
correspond to observed equilibrium shapes at temperature T = only; however, because the energy scale Kb is large 
on the scale of fc^T, thermal fluctuations are generally negligible and play an important role only for "soft" modes, 
especially near instabilities (Wortis et al., 1997). The Euler-Lagrange equations are in general nonlinear and can 
have multiple solutions. The lowest-energy solution for given A and V is automatically stable to small fluctuations; 
higher-energy solutions must be tested for stability. All local-minimum solutions are candidates for stable observable 
shapes, except in exceptional cases (near instability) where energy barriers become comparable to ksT. 

The stable shapes produced by this process depend on the control parameters: the geometrical parameters, A and 
V; the determinants of curvature, Cq and A^qj which only enter in the combination (0) or, equivalently, (|^)-(^); the 
moduli, K, K, Ka, and /i; plus, finally, any parameters required to characterize the relaxed shape of the membrane 
skeleton. The area and volume of a typical RBC we take as (Bessis, 1973) A = 140^m^ and V = QO/im"^. The moduli, 
given in Sec. I, are less well determined by experiment (we shall have more to say on this point in Sec. V). This 
leaves as unknown control parameters the curvature variables and the cytoskeletal parameters. 

The Euler-Lagrange variational equations derived from ( |l^ are not numerically tractable except in the case of 
axisymmetry, which clearly does not apply to echinocytic shapes. In this paper, we aim to treat only fully developed 
echinocytic shapes (echinocyte III). For this case, individual spicules are identical and axisymmetric in shape to a 
good approximation and, in addition, the central body which they decorate is approximately spherical with radius 
Rq. The observed distribution of spicules is rather regular, and we shall approximate the local spicule packing as 
a triangular array (except for special numbers of spicules, there must, of course, be some defects), which will look 
increasingly like an hexagonal close-packed structure, as the number of spicules, n^, becomes large. 

The base of each individual spicule, where it meets the sphere tangentially, is a circular contour L of radius r^. If 
9l is the angle subtended by L at the center of the sphere, then 

rL^RosinOL- (13) 

Because of the close-packed structure, the circles L meet tangentially. We explain in Sec. Ill how to derive appropriate 
boundary conditions where the spicules meet the sphere. Solving the axisymmetric Euler-Lagrange equations then 
determines a spicule shape, including an individual spicule volume Vs and area Ag. In terms of these variables, the 
overall area and volume of the RBC are taken to be 



A « nAs + 0.1 X AttRI (14) 

and 

Air 

V^—Rl+nVs, (15) 

with the spicule number 

n, «3.6 X (^y^^ . (16) 

In writing these relations, we have assumed that 10% of the spherical surface is not covered by the circular spicule bases 
(close packing on a flat surface would give 9.3%; curvature effects and packing defects both increase this number). 
The spicule volume Vg is calculated with respect to a plane through L, and Eq. |l5| neglects curvature corrections. 
These approximations are good when the number of spicules is large, as it will turn out to be for the echinocyte 
shapes. 
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III. THEORY 



A. Membrane mechanics 



In our treatment, spicules are assumed to be axisymmetric, and individual spicules are joined to the central body 
along the contour L, as illustrated in Fig. 2. Thus, our calculation involves finding a family of energy-minimizing 
axisymmetric spicule shapes and selecting from that family those shapes consistent with appropriate mechanical 
boundary conditions applied along L. 

Parametrization of the axisymmetric spicule shape is illustrated in Fig. 2. The variables z and r measure distances 
along and perpendicular to the symmetry axis, respectively, while s measures the arclength from the pole. The 
function z(r) determines the shape. 9 is the angle between the local normal and the symmetry axis; Cm{r), Cp{r) are 
the principal curvatures, 

^ d& , ^ sin 61 , „, 

C„ = — and Cp = . 17 

as r 

In calculating the spicule shape, we shall assume that the relaxed cytoskeleton is locally flat, a good approximation as 
long as the number of spicules is large, so that the spicule size is small compared to Rq. Thus, in forming each spicule, 
a flat circular patch of relaxed cytoskeleton must be elastically deformed to fit the spicule contour. We assume that 
this deformation is axisymmetric, so that the center of the patch remains at the apex of the spicule, and each point 
of the patch at relaxed radius Sq maps to a point at arclength s and radius r on the spicule contour. The principal 
extension ratios can then be written, 

Ai = -, X2 = ^. (18) 
So aso 

When these expressions for Ai and A2 are inserted into Eq. ^ we obtain an explicit form for the elastic energy in terms 
of the functions so{s) and r(s), 

f IT ds 
Fci = nKa / sodso ; 1 



So dso 

+n,fs,dso(^'-^+'-^^~2). (19) 
J \so ds r dso J 

Using Eqs. |, |, |l^, and |l|, we can implement the stationarity condition for the free-energy functional ( p^ ) with 
respect to variations of membrane shape and cytoskeletal strain. This leads to a set of five coupled first-order 
differential equations, which are written down explicitly in the Appendix. Note that the ADE term in Eq. ^ enters 
only through the appearance of an effective spontaneous curvature 

Cf^Co~^{AA~AAo), (20) 

which must be determined self-consistently via Eq. [sj Although these Euler-Lagrange equations look complicated, they 
can be related to simple mechanical force-balance conditions in a way that makes their content entirely transparent 
(Mukhopadhyay and Wortis, in preparation). We have used these equations to calculate (approximate) spicule shapes. 

We now turn to the boundary condition where the spicules meet the sphere (see Fig. 2). Along each element of L, the 
spicule membrane exerts a tension r[| in the plane of the membrane and perpendicular to L, and a tension t± normal 
to the plane of the membrane. There is no tension in the third ("shear") direction because of the membrane fluidity. 
In addition, the membrane skeleton exerts an independent tension which must be in-plane, since the skeleton lacks 
bending rigidity. In general, these tensions would be different at different points of L; however, in the axisymmetric 
approximation, the tensions are uniform around L and the skeletal tension is directed radially. To maintain mechanical 
equilibrium, these tensions must balance where the spicules meet at point A. Note that A is a point of bilateral 
symmetry between the two adjacent spicules, so that in-plane tensions, which act in opposite directions, always 
balance by symmetry. On the other hand, the normal tensions t± must vanish individually, since they act in the same 
direction. The normal tension is related to the isotropic bending moment per unit length, 

M=Kb{Cp + an-Cf), (21) 

by (Evans and Skalak, 1980; Mukhopadhyay and Wortis, in preparation) t± ~ dM /ds. Thus, an appropriate boundary 
condition along L is 
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"^^^"ds '^"'^ " ^^^^ 
This is the boundary condition which we shall apply in our calculations. Of course, this boundary condition is 
approximate, since the real spicule is only approximately axisymmetric. Indeed, with equal logic, we could argue by 
symmetry that t± should vanish at point B, where three adjacent spicules meet at a radius from the spicule axis some 
15% larger than r^. We shall use this observation in Sec. IV to get a rough measure of the error introduced by the 
approximation (|2^). 

We now outline the algorithm used to solve numerically the Euler-Lagrange equations for stationary spiculated 
shapes. A "shooting method" was employed to integrate the Euler-Lagrange equations from the apex to the edge 
of the spicule, the contour L along which the normal tension t± vanishes. The set of solutions can be characterized 
by five parameters. Three of them are the global parameters: the pressure difference P, the Lagrange multiplier E, 
and the effective spontaneous curvature Cg*^, Eq. ^ The remaining two parameters correspond to initial values 
of variables in the integration procedure, i.e., the curvature (Cm — Cp) and the stretching ratios (Ai = A2) of the 
cytoskeleton at the pole. Integrating the Euler-Lagrange equations to a point where Eq. 22 is satisfied determines 
Tl and 6^ and, therefore, Rq and Ug from Eqs. |l^ andjl^. The initial conditions may be adjusted iteratively to fit 
four parameters, A (Eq. |l^), V (Eq. 15), AAq (Eq. and the cytoskeletal stretching ratio, 7, which measures the 
ratio of the area of the calculated spicule to that of the corresponding relaxed cytoskeletal patch (roughly the same as 
the ratio of the area of the full echinocyte to that of the full relaxed membrane skeleton, since the corrections in Eq. 
[1^ for the spherical segments are small). We are left, finally, with a one-parameter family of mechanical-equilibrium 
solutions, which we take to be labelled by the spicule number n^. The predicted equilibrium configuration is the 
solution nf^J^ which minimizes the overall energy. Of course, for the exact solution, Us is a discrete variable: Probably, 
there exists, for given initial conditions, a set of distinct branches of solutions characterized by different rig's and 
different arrangements of spicules on the surface, as discussed more fully in Sec. V. 



B. Scaling analysis 

Before proceeding to solutions, we wish to identify the important energy and length scales of the problem. Because 
we are working at a temperature which is effectively low, one overall energy scale drops out of the mechanical problem. 
We take this scale to be set by the bending modulus (nh and k are similar in magnitude). There are three length 
scales: The first is the overall scale of the RBC, which we denote R. The second, we may take to be I/Cq from 
Eq. (which is equivalent to using a length based on AAq, Eq. ||); alternatively, we can take the second length to be 
1/Cq, where Cq^ is the effective spontaneous curvature from Eq. EO (we assume that Co and Cg* are comparable 
in magnitude). The third length is the elastic length scale A^i, Eq.lT. In general, there may be other length scales 
associated with the relaxed cytoskeletal shape; however, we shall assume that any such lengths are comparable to R. 
For the RBC, R ^ A^i- The length scale I/Cq is a control parameter whose magnitude can be tuned across the range 
defined by R and A^i, resulting in the observed RBC shape classes. 

The size of the RBC sets the total area and volume in the problem. Note that the RBC has a relatively smaller 
volume for a given surface area than a sphere. If a sphere with the same surface area has a volume l^phere, then the 
ratio of RBC volume to that of the equivalent sphere defines a reduced volume, v = Vrbc /Vsphcrc- For the RBC, 
V w 0.6, significantly less than unity, which allows the RBC to assume its sequence of distinct nonspherical shapes. In 
its usual discocytic phase, the red blood cell is expected to have a magnitude of 1/Cq that is of the order of the red 
blood cell dimensions or larger. Tuning Cq to sufficiently large positive values generates the sequence of evaginated 
echinocytic shapes; tuning it to sufficiently large negative values generates the sequence of invaginated stomatocytic 
shapes. 

The significance of the elastic length scale A^; is that it determines how strong the effect of elasticity is for structures 
with some characteristic length, say L. To see this, we imagine rescaling the problem with the energy scale ki, and the 
length scale L in such a way as to make all the parameters appearing in the Euler-Lagrange equations dimensionless. 
In this way, we arrive at the rescaled parameters, 

Kf, = 1 

P' = PL^/Kb 

H = fJ-L'^/nb 

Cf = CfL. (23) 



8 



Equivalent rescaling of the variables, C„j = CmL, etc., leaves the form of the Euler-Lagrange equations unchanged, 
although the parameter values (Eq. p3|) vary with L. For example, suppose L is the echinocyte spicule radius. If 
L > Aei, then /i is greater than unity, and the elastic terms from Eq. ^ dominate the bending-energy terms; conversely, 
if L < Ae;, then /i is smaller than one, and the bending-energy terms dominate. Thus, generically, elasticity is strong 
at length scales larger than A^i and weak at length scales smaller than Ag;. This statement may seem surprising 
to those who think of the discocyte as effectively shaped by bending energy. The key point here is the plasticity 
of the membrane skeleton, albeit on time scales much longer than the induced shape changes we are discussing: If 
the relaxed membrane skeleton closely matches the shape preferred by bending energy alone, then (of course) elastic 
effects are small. This is apparently the case for the resting discocyte (but definitely not for the echinocyte!). 

A similar scaling argument gives insight into the characteristic spicule size in the echinocyte region. Suppose we let 
L — Aei in Eq. At this length scale, the bending and elastic terms (k^^, K^, and fi ) are comparable and the only 
remaining parameters are P , S , and Cg^ (or, equivalently, Cq). If, in addition, P and E are small, then the Euler- 
Lagrange equations are equivalent to those for the unconstrained minimization of the combined bending-plus-elastic 
problem (H). But, this is a problem which we understand: When Cg^ is small in magnitude, then the membrane 
remains flat (at the length scale of A^i); on the other hand, when Cq^ is large in magnitude, then bending energy 
(which dominates the shape at scales below A^i) will favor budded shapes on the scale of l/Cg*^, either invaginated 
or evaginated depending on the sign of Cq^ . Between these limits, when Ae; and are comparable, we expect 

spicules (for positive Cg*) scaled by tending towards smoother shapes when Aei is somewhat smaller than 

and towards sharper, more columnar shapes when Aei is somewhat larger than 1/Cq^. Note that, as long as 
we rescale to the length A^i, the expected spicule shape is predicted to be independent of cell size R and to depend 
only on the scaled value of Cg^ and the ratio of the two elastic constants, n and Ka- 

The above argument depends on the condition that the scaled values of P and S be small at the elastic length 
scale. For smooth and flaccid shapes like the discocyte (and when the cytoskeleton is not significantly stretched or 
compressed) the pressure difference P is generated by the bending energy, so P is of order unity on the scale of the 
cell size R, and similarly for E . It then follows from Eq. ^ that the condition is well satisfied for the discocyte and 
other nearby shapes (as we shall see, the initial echinocyte shapes fall into this class). Of course, if the system is 
pushed too hard, then this condition breaks down and other length scales enter the spicule-shape problem. 

In summary, the shape of a spicule arises as a result of the interplay between the two length scales. Ad and . 
When is positive and much larger than Ad, the effect of elasticity is strong and we can expect spicules that are 

at most gentle bumps. It is only when becomes of the order Aei or smaller, that we may expect sharp spicules. 

Thus, if we look at the series of shape transformations, from discocytes through discocytes with gentle bumps to 
sharply spiculated structures, the sharp spicules will first arise with a radius comparable to Aei. We shall verify this 
scenario in the next section by explicit solution for the stable shapes. 

IV. RESULTS 

In this section we report echinocyte shapes calculated according to the program outlined in Sec. Ill and using 
the elastic moduli (/«{,, k, /i, and Ka) given in Sec. I and the RBC parameters {A and V) given in Sec. II. The 
cytoskeletal stretching ratio 7 is not known to good accuracy for the RBC. In this section, we take 7 = 1; however, 
we have tested values in the range 0.7-1.2, and it is one of our results that variation within this range produces no 
appreciable changes in the number of spicules or their shape. The spontaneous curvature and relaxed area difference 
are not known and presumably vary somewhat over a typical RBC population; however, they enter the mechanics 
problem only in the combination Aoq, Eq. |9[ which we take as our principal control parameter. 

Figure 3 displays the equilibrium spicule shapes that we obtain for a sequence of values of Aoq. Figures 4 and 5 plot 
the corresponding calculated number of spicules which appear on the fully developed echinocyte III at equilibrium 
and the corresponding values of Cq^ as functions of Aag. 

How well do the calculated shapes and spicule dimensions agree with experimental observations (Bessis, 1973; 
Brecher and Bessis, 1972; Chailley et al., 1973). As shown in Fig. 3, the calculated shapes fall into three distinct 
categories. For Ug less than about 20, corresponding to Aoo appreciably smaller than 0.016, the spicule shapes are 
broad and rather hat-shaped, rather different from those seen in experiment. As Aqq and Cq^ increase, the number 
of spicules increases and the individual spicules become smaller, just as predicted by the scaling argument at the end 
of Sec. III. By the time Aao is at or above 0.016, for Us in the range 30-60, the spicule shapes are columnar with 
rounded tops, in good general agreement with the kinds of shapes seen in experiment. Beyond Ug = 70, the spicules 
become narrow-necked, a shape not typically seen in experiment but illustrating the dominance of the bending energy, 
as discussed in Sec. III. 



9 



In the "good," central region, the spicule dimensions and numbers are very comparable to those seen in experiment. 
A typical example (Aao = 0.018) is shown in Fig. 6(a). For this case, we find = 41, a radius of the central, 
spherical body of Rq = 2.57, and an aspect ratio, height over width at base (which also corresponds approximately 
to the distance between adjacent spicules), of about 0.8, all of which are in reasonable agreement with observation. 
The only significant discrepancy is that the spicule height is close to 1.0 /um, somewhat larger than the range 0.6-0.8 
/xm seen in experiment, as is the aspect ratio. We shall comment on possible reasons for this discrepancy in Sec. V. 
We emphasize that these shapes are a prediction of the model, not an assumption as they are in other recent work 
(Iglic, 1997; Iglic et al., 1998a, b). The only other spicule-shape predictions are those given by Waugh (1996), whose 
calculation suggests cross sections that arc less steep at the sides and more pointed at the apex than what is observed. 
Whether this is a consequence of the variational form assumed or of the neglect of stretching elasticity in the model is 
not clear. The calculations of Iglic et al., based on a postulated spicule shape and assuming skeletal incompressibility 
{Ka oo), produce somewhat higher values of the aspect ratio. At a given reduced area difference Aao, they find a 
spicule number which is larger than ours by almost a factor of two. We will return to some of these issues in Sec. V. 

Before discussing further what we shall argue arc the non-biological regions of hat-like and narrow-necked shapes, 
we must review the full echinocyte sequence (Brecher and Bessis, 1972; Bessis, 1973; Chailley et al., 1973; Mohandas 
and Feo, 1975). Experimental observations of RBC's which are initially smooth, axisymmetric discocytes show that 
echinocytosis occurs in three stages. The first stage (echinocyte I) corresponds to irregular disc-like shapes with a 
central dimple. In the second stage (echinocyte II) the irregularities evolve into more or less well defined spicules 
and the central dimple disappears, resulting in an oblately elliptical body decorated with spicules. In the third stage 
(echinocyte III) the central body becomes approximately spherical and the spicules become somewhat smaller and 
more numerous. When pushed beyond the echinocyte III, some kind of disconnection between the plasma membrane 
and the membrane skeleton occurs (see below) and plasma membrane is shed from the spicule tips into microscopic 
exovesicles. The remaining spherical shape with pointed protuberances is called a sphcrocchinocyte (Bessis, 1973). 

Our calculation explicitly assumes that the central body is spherical and is, therefore, appropriate only for the 

late stage-Ill echinocyte. What we believe happens is that, for Aao ~ 0.016, the branch of echinocyte III shapes 

becomes energetically unstable to a branch of echinocyte II shapes. There is no internal way that our calculation 
can show that this crossover takes place; however, there are two pieces of evidence we find convincing. The first is 
a simple estimate: A separate calculation of the discocyte branch shows that its shape (and AA, for example) is 
roughly independent of Aao, allowing us to estimate its energy as a fimction of Aao. The calculated energy of our 
echinocyte III branch is lower than the discocyte branch but increases as Aao decreases. The crossing point of these 
two branches at Aao ~ 0.016 is, thus, an approximate lower limit to the transition out of the echinocyte III phase. 
The actual transition out of the echinocyte III branch must presumably take place at a somewhat higher value of 
Aao. The second piece of evidence is that, in separate work which we shall report elsewhere (Lim, Wortis, Boal, and 
Mukhopadhyay, 2002), we have done numerical simulations of RBC shapes using the same set of numerical parameters 

chosen here, and we find just such a transition at Aoq ~ 0.017 from the echinocyte III to a branch of echinocyte II 
shapes. 

At the opposite end of the range of "good" spicule shapes, our calculation fails again and, again, for reasons that 

cannot be assessed internally. Here the size of the predicted spicule-set by the scale of 1/ Co*- becomes comparable to 
the length of the elementary spectrin tetramer of the membrane skeleton (Liu et al., 1987). At this length scale, the 
continuum picture breaks down for the membrane skeleton (although it remains valid for the lipid component). What 
presumably happens is that a point of instability is reached, where lipid flow occurs and the plasma membrane buds 
outwards in regions between the cytoskeletal anchors. Such buds, lacking cytoskeletal support, are known to be fragile. 
They presumably break off. leading to microvesiculation, membrane loss, and the observed sphcrocchinocytosis. 
Naturally, this vesiculation is expected to happen first at the spicule tips, where the elastic network is maximally 
expanded and the distance between anchors is largest. 

In nature, the plasma membrane and its cytoskeleton are bound together by protein anchors; in our model this 
binding is replaced by the condition that the bilayer and the elastic network co-exist on the same mathematical 
surface. This requirement means that there is a local force per unit area or pressure Q with which the skeleton pulls 
inward (Q > 0) or pushes outward {Q < 0) on the membrane surface (or, equivalently, with which the membrane 
pulls outward or pushes inward, respectively, on the skeleton). The expression for this pressure is (Mukhopadhyay 
and Wortis, in preparation) 
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and it is plotted in Fig. 7 as a function of the radial distance r from the spicule axis for Aoq = 0.018. As expected, Q 
is positive close to the tip of the spicule and negative around the neck where the skeleton has been compressed. The 
magnitude of the pressure at the spicule tip is of the order of 20 pN//xm-^. It is useful to compare this pressure to the 
critical disjoining pressure necessary to break the anchoring and to separate the membrane skeleton from the bilayer. 
This critical pressure has not, to the best of our knowledge, been measured; nevertheless, it can be estimated. We 
take 10 pN as a crude estimate of the force required to extract quickly a single anchoring protein from the bilayer. 
Taking into account the density of anchors, we estimate a detachment pressure of about 2000 pN///m^ (consistent 
with measurements reported in Waugh and Bausserman, 1995), two orders of magnitude larger than what we have 
calculated above. This comparison is important, as membrane-cytoskeleton disjoining is another potential mechanism 
for spheroechinocytosis, an alternative to the one discussed in the paragraph above and has, indeed, been proposed by 
Iglic and others (Iglic et al., 1995, 1996). We conclude that such direct disjoining seems unlikely, unless it is associated 
with anchoring defects. 

We have seen for the human erythrocyte that the fully developed echinocyte III, with well developed spicules, 
first appears when l/Cg*^ becomes small enough to be comparable to the elastic length scale Aej. It is a plausible 
hypothesis that this is a general connection. If so, there is a direct relation between the spicule dimension at the 
onset of the echinocyte III and the ratio (Eq. |l|) of the bending rigidity to the elastic moduli. Interestingly, Smith et 
al. (1982) have studied spiculated shapes for the red cells of a variety of species. Despite considerable variation in 
the volume of the red cells-ranging from the goat RBC, which is appreciably smaller than the human one, through 
the elephant seal, which is appreciably larger~the spicule size is quite stable, suggesting that the elastic length scale 
is not strongly variable from one species to another. This means, of course, that the number of spicules on the typical 
echinocyte is small for the goat 10) and large for the elephant seal (~ 100), as observed. To illustrate this point 
further, we plot in Fig. 8 spicule shapes at different volumes and different ratios Ka/n but with Ag; held fixed. 
We find that at smaller volume, the spicules sizes do become smaller but the dependence on volume is weak. The 
dependence on the ratio of elastic constants K^/ n is even weaker. At a value of that is larger by a factor of 2 
{Ka = 6/x), the spicules become somewhat broader and their number decreases by around five percent. Of course, the 
local stretching of the cytoskeleton both at the spicule apex and in the neck region close to the base may be expected 
to change considerably. 

These results are based on the model described in Sec. Ill, which assumes spicule axisymmetry, as incorporated in 
the (approximate) boundary condition, Eq. In closing this section, it may be useful to comment on the reliability 
of this central assumption. We have two remarks. The first concerns the boundary condition. The neighborhood 
of a given spicule is clearly not axially symmetric, so this condition is clearly approximate. If the calculated shapes 
depended sensitively on the precise way that this boundary condition is applied, then results based upon it would be 
suspect. We motivated the boundary condition {t± — 0) by looking at the point labelled A in Fig. 2, where adjacent 
spicules are tangent. With equal logic, we could have argued this boundary condition at the point labelled B. If we 
made this choice, how different would the calculated spicule shapes turn out to be? We have run some test cases 
and find changes in, e.g., the spicule height of less than 1%, suggesting that the results are strongly insensitive to the 
details of the boundary condition. But, this is far from sufficient. As a second test, we have recently completed a 
program (Lim et al., in preparation) of direct simulation of the energy functional (j2|) on a triangulated surface. This 
program allows us to calculate echinocyte shapes without assumptions about spicule axisymmetry. We find echinocyte 
shapes in this region of parameter space and can compare the simulated spicule shapes with those calculated here at 
the same parameter values. Figure 6 shows one such comparison. Although the base of the simulated spicule is a little 
narrower than that given by the approximate axisymmetric calculation, the overall level of agreement is excellent. At 
these parameter values, the number of spicules in our calculation is 41, while the number found by simulation is 34. 
Similarly, the radius of the central spherical body to which our spicules are fitted is 2.57 /xm, while that found in the 
simulations is about 2.5 /im (the central body is only roughly spherical in this case). The approximate calculation 
presented in this paper appears to be reasonably reliable. 

V. DISCUSSION 

The calculations presented here build on the previous work of Iglic et al. (Iglic, 1997; Iglic et al., 1998a, b) and show 
that a simple mechanical model of a uniform composite membrane can give a good account of observed echinocyte 
shapes. The central control parameter in this process is the effective area difference (or curvature) of the lipid 
component, in agreement with the bilayer-couple hypothesis of Sheetz and Singer (1974). The role of the cytoskeletal 
elasticity is to suppress the formation of narrow-necked buds which would form in the absence of cytoskeletal shear 
resistance. Does this good agreement prove that mechanics (as opposed to biochemistry) is the sole determinant of 
echinocyte shape? Of course, not; however, it does strongly suggest that simple mechanics will play an important role 
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in any complete picture. We end this section with a brief discussion of some loose ends and future directions. 

A. Elastic constants and elastic nonlinearities 

There has been much discussion in the literature of what are the "correct" elastic constants to use for the membrane 
skeleton. We have chosen /i « 2.5 x 10~^ J/m^ and ~ 3/i. These values are in rough agreement with recent 
measurements at low deformation (Henon et al., 1999; Lenormand et al., 2001) and somewhat lower than those 
reported in connection with pipette aspiration of red cells (Evans, 1973a, b; Waugh and Evans, 1979), where large 
deformations play an important role. 

Discussion of values of the elastic moduli is complicated by the occurrence of nonlinearities in the elastic response. 
At low deformation. A; = 1 + with ti small, there are only two rotational invariants from which to construct the 
elastic energy, \{ei +£2)^ and —£2)^, whose coefficients define the linear moduli of stretching and shear elasticity. 
It is easy to verify that these linear moduli are just Ka and /i, respectively. However, at large deformation, terms of 
higher orders in (ei + £2) and (£1 — £2)^ will generally also appear in the elastic energy, each with its own coefficient. 
The expression (H) for the elastic energy makes a particular choice for each of these coefHcients, one which is simply 
proportional to the linear moduli. At small extension ratios, these higher-order terms are unimportant; however, 
for the echinocyte shapes, we find ti of order unity, so these terms cannot be neglected. Unfortunately, present 
experiments do not constrain the values of the higher order coefficients significantly, although it is probably not an 
accident that the moduli derived from the pipette aspiration experiments (Evans, 1973a, b), which have extension 
ratios up to nearly 2 (Lee et al., 1999), are larger than those measured in the linear regime (Henon et al., 1999; 
Lenormand et al., 2001). This suggests a hardening of the elasticity at large deformation, consistent, for example, 
with terms in (£1 + £2)'^, (£1 + £2)(£i — £2)^, and so forth, coming into play. We will deal with the issue of nonlinearity 
at more length in a future publication. For the present, we simply remark that realistic spicule formation does require 
a proper balance between stretching and shear energies at large deformation. Our chosen ratio of Ka/ ^ — 3 suffices; 
however, a lower ratio plus some additional nonlinear hardening would also achieve the same effect and may be more 
realistic. Any effect that changes spicule shape significantly may also be expected to modify quantitatively the regions 
of stability of different shape classes. 

Previous variational work on echinocyte shapes by Iglic and coworkers (Iglic, 1997; Iglic et al., 1998a, b) has 
assumed local incompressibility of the cytoskeleton {Ka — > 00), thus forcing all the elastic energy into the shear term. 
This was done for practical reasons; we are now in a position to understand the effect of this approximation. Figure 
8(c) suggests that increasing Ka at fixed ^ has a modest effect on spicule shape, decreasing the height and the base 
radius somewhat without changing the mean radius appreciably. This is consistent with the generally larger number 
of spicules predicted by Iglic et al. (Iglic, 1997; Igfic et al., 1998a, b). Because of the fimited shapes available in the 
variational parametrization, these authors cannot follow the full spicule-shape evolution exhibited in our Fig. 3. 

B. Spicule placement, metastability, and related matters 

Our calculation assumes for the echinocyte regular spicule placement with six-fold coordination on a central spherical 
body. Of course, this is only an approximation. First, there is the topological requirement for a net undercoordination 
of 12. More generally, if the energy minimization problem were to be solved exactly, we would expect to find several 
distinct sheets of locally-minimizing shapes with different spicule numbers Us and different spicule organization. For 
any given set of parameter values, one such shape would have the lowest energy; however, several others might have 
nearby energies. Metastability of a shape would be expected, since the energy landscape is on the scale of k;,, which 
is large on the scale of kBTroom- This metastability, in turn, would cause hysteresis in the shape transformations. 

Some experiments show that, as the red cell is chemically cycled back and forth between discocyte and echinocyte 
phases, spicules always appear at the same locations on the cell surface (Furchgott, 1940; Bessis and Prenant, 1972), 
suggesting that spicule placement could be intrinsically related to defects and other inhomogeneities in the cytoskele- 
ton. This would not happen in a continuum model. Without wishing to dispute the experimental evidence, we 
can only say that our work shows that spicule formation does not require cytoskeletal inhomogeneity. Of course, if 
inhomogeneities exist, then they will have an effect in selecting the pattern of spicule placement. 

It is not clear a priori whether or not the spicule is a solitonic object with its own intrinsic scale and, if so, whether 
spicules tend to attract or repel at long distance. Some recent work (Lim et al., in preparation) has suggested long- 
range attraction, which would mean that, under appropriate circumstances, spicule clumping might occur, leaving 
"bald" regions on the echinocyte surface, as is seen in some pictures of echinocyte shapes (see, for example, Iglic et 
al., 1998a). None of this finer detail is captured in our approximation. 



12 



C. Extended shape/phase diagram 



We have focussed on the spontaneous curvature or area difFerence Aag as the principal parameter controUing the 
evolution of echinocyte shape. Of course, induced changes in any parameter of the membrane energy expression Eqs. 

and the constraints V and A, singly or in combination, will affect the cell shape. A (bio)chemical modification 
of the red-cell environment presumably changes all these parameters to some extent; it is just that some changes 
are larger and more important than others. If there are situations in which laboratory reagents cause important 
modification in other parameters-e.g., the elastic constants-then addition of these reagents will cause the red cell to 
be driven along a trajectory in an extended parameter space. In such a situation, we would need to consider the 
minimum-energy shape (or shapes) as functions of several variables in an extended shape/phase diagram. The work 
surrounding Fig. 8 illustrates the beginning of such a study. 

We have summarized in Sec. I the mechanism by which some common reagents affect the Sheetz-Singer (1974) 
parameter Auq. The pH effect and the associated glass effect are less readily related to Aao but may well be primarily 
related to another dimension of the generalized phase diagram. For example, Gedde et al. (1995) have shown that 
lipid asymmetry and the presence of inner leaflet titratable groups, which might be expected to be associated with 
Aoo, do not have appreciable influence on the pH effect; however, Elgsaeter et al. (1986) and Stokke et al. (1986) have 
shown that the membrane skeleton expands in vitro in response to high cytoplasmic pH. Now, according to the scaling 
argument (Eq. pi] ), we know that such expansion is equivalent to decreasing the shear modulus /i and, therefore, to 
increasing the elastic length scale A^i (Q). But, echinocytosis occurs when l/Cg^ shrinks below A^i. Thus, at fixed 
Cq* (roughly equivalent to fixed Aao), expanding the membrane skeleton could promote echinocytosis, as observed. 
It would be premature to claim that this is the full explanation of the pH effect. The point is that, as long as changes 
in biochemical variables produce homogeneous changes in the material parameters that appear in Eqs. our model 
continues to apply, although the relevant experimental trajectories may be in a larger parameter space. 

D. Membrane composition, inhomogeneity and charge effects 

Description of any effects which lead to inhomogeneity in membrane properties requires modification of our model. 
One such effect which is certainly present arises from the inhomogeneous distribution of lipids and/or proteins. We 
know that the plasma membrane is composed of a mixture of lipids and proteins. The bending modulus kj, depends 
on membrane composition. As long as the composition remains homogeneous, the model Eqs. ^ ^ holds; but, 
any mechanism that produces compositional inhomogeneity on the scale of RBC shape features would necessitate 
modification of the model. Thus, for example, observed raft formation (Simons and Toomre, 2000; Pralle et al., 2000) 
suggests membrane inhomogeneity on the scale of 0.05-0.07 /xm. When local radii of curvature become comparable to 
the inhomogeneity scale, then inhomogeneity is expected to influence shape. This may occur for the echinocyte near 
the spheroechinocyte limit; however, there is no evidence that it is a dominant effect elsewhere. 

Other interesting effects involve the coupling of lipid composition to geometry. Thus, for example, outer-leaflet lipids 
with large heads may be expected to segregate preferentially to spicule tips, as would inner-leaflet lipids with large tails. 
Such an effect could be included in a mechanical model by adding a composition field (or fields) and including terms 
coupling composition to curvature. Alternatively, one might imagine mechanisms involving membrane components 
which preferentially favor (or avoid) regions around the cytoskeletal anchoring complexes. These components would 
be relatively dilute (concentrated) in regions where the cytoskeleton was significantly expanded (as it is near the 
spicule tip). To represent this effect, one would need to couple composition to local cytoskeletal strain. 

While these mechanisms may be present to some extent in the red cell, our calculations suggest they are not in any 
way required for spicule formation. 

Finally, chemical groups which carry net charge at physiological pH are common both for lipid headgroups and 
for cytoskeletal components. As long as these charges are compensated by counterions on scales smaller than the 
local radii of curvature, they enter the mechanics only via their effect on the mechanical moduli. On the other hand, 
when relevant lengths become comparable to or smaller than the Debye length, then charge effects must be treated 
as nonlocal. 



E. Conclusions 

Our calculation provides echinocyte sizes and shapes in excellent agreement with experiment. Although shapes are 
calculated under the approximation of spicule axisymmetry, there is evidence (e.g., Fig. 6) that this approximation 
is good. The approximation does not allow us to estimate internally the range of echinocyte III stability; however. 
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extrinsic arguments (Sec. IV) suggest that observable shapes should be bounded at low Aaq in the region where ^/C^^ 
first becomes smaller than Agj and at large Aao in the region where the spicule dimension becomes comparable to 
the cytoskeletal discreteness. And, indeed, outside this region the predicted shapes are no longer seen experimentally. 
While this consistency lends credibility to the continuum picture, the way in which biochemical influences drive 
the control parameters of the model is not addressed. The Sheetz-Singer parameter, Aao, which measures the area 
difference between the lipid monolayers (and the related effect of membrane spontaneous curvature) does seem to be 
the dominant control parameter for some kinds of induced shape change. On the other hand, other effects-such as the 
pH effect-may probe other dimensions of the shape/phase diagram, e.g., the role of the cytoskeletal elastic constants 
and prestress. We hope that having a mechanical model with predictive power may eventually help to elucidate 
important biochemical questions. From this perspective, experiments can now focus on the way in which biochemical 
probes affect all the mechanical control parameters, including elastic and cytoskeletal variables, in addition to Aaq. 
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VI. APPENDIX 

The Euler-Lagrange equations obtained by minimizing $ may be written down as the following set of five coupled 
ordinary differential equations expressing the variables sq, r, 6, Cm and h as functions of the arclength s: 

f ds r \ j So r 
2K^r 1 + /iso — 




ill 

+Hbr [Cl - 2CfCp - Cl,] - h{s) cos 9 + Pr^ sin 6* = 0, (25) 

^=cose (26) 
ds 

^ = Cm (27) 



ds 

dCm cos 9 sin 9 C,r, cos 9 Pr cos 9 b sin 9 



ds r r 2Kb 2rKb 

db f ds r 

2K^[- 1 



(28) 



ds \dso So 



dsQ f 1 dsQ s ds . 
-~^^so— 3 5-3— + 2a 



ds \ Sq ds dso ^ 

+Kb [cl, - 2CfCm - Cl] + 2pr sin 9. (29) 

Here b{s) is a local variable introduced to impose the constraint, dr/ds — cos 6*, as first suggested by Peterson (1985). 
The variables s, sqj f^ and b vanish at the North pole. Note also that the surface tension a used in these equations 
is, in general, shifted from the parameter E appearing in Eq. 
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FIGURE CAPTIONS 



Figure 1 (a) A budded or vesiculated shape, which would be a low-energy configuration of a lipid membrane with- 
out cytoskeleton at sufficiently positive spontaneous curvature, (b) The spiculated shape into which (a) is 
transformed by the elastic-energy cost of cytoskeletal deformations in the high-shear neck regions. 

Figure 2 Parametrization of a single-spicule shape, showing the boundary L where the spicule joins the spherical 
central body in constructing the full echinocyte shape. The bases of two adjacent spicules are sketched in 
to illustrate the points A of spicule tangency and the point B which is symmetrically situated between three 
neighboring spicules. 



Figure 3 Sequence of spicule shapes for increasing values of the reduced effective relaxed area difference Aao, Eq. 




which measures the combined effect of spontaneous curvature and additional area in the outer leaflet in driving 
the membrane to bend outward, (a)-(d) correspond to Aao = 0.014, 0.018, 0.020, and 0.022, respectively. Note 
the bar shows the elastic length scale Aej. Note how the spicule sharpens as Aao increases and how the neck 
begins to form when the spicule dimension (set by 1/Cq*) falls below K^i so that the bending energy begins to 
dominate. 

Figure 4 Calculated equilibrium number of spicules Ug plotted as a function of the reduced effective relaxed area 
difference Aao. 

Figure 5 Calculated effective spontaneous curvature Cq*, Eq. (po[), plotted as a function of Aao. 

Figure 6 Calculated spicules at Aao = 0.018. (a) Shape found in the axisymmetric approximation, as discussed in 
this paper (n^ = 41, i?o = 2.57/im). (b) Shape found by numerical minimization on a triangulated surface 
{us = 34, i?o = 2.5/im). 

Figure 7 The local pressure Q, Eq. (p^), exerted by the cytoskeleton on the bilayer plotted as a function of radial 
distance, for a spicule corresponding to Aao — 0.020. 

Figure 8 Dependence of spicule shape on parameters at Aao = 0.018. (a) Shape with the standard RBC parameters, 
as discussed in the text, (b) Effect of starting from (a) and increasing by a factor of 2 (i.e. — 6/i). 
The spicules become slightly broader, (c) Effect of reducing overall cell dimension by a factor of two, so that 
V — > V/% and A — > A/4,. This change reduces spicule size by about 10%. 
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FIG. 3. 
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FIG. 4. 
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FIG. 5. 
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FIG. 6. 



25 







■10 

2 4 6 8 

Radial distance (x 0.1 microns) 



FIG. 7. 



26 




FIG. 8. 
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